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THE  DEVELOPMENT  OF  A  DYNAMIC  ELASTOPLASTIC 
FINITE  ELEMENT  ANALYSIS  FOR  FAST  FRACTURE 
UNDER  IMPACT  LOADING 


by 


Jalees  Ahmad  ,  J.  Jung,  C.R.  Barnes,  and  M.F.  Kanninen 


INTRODUCTION 


Fracture  mechanics  researchers  are  becoming  aware  that  the 

applicability  of  even  rigorous  dynamic  analyses  of  unstable  crack  propagation 

and  crack  arrest  may  be  more  limited  than  was  previously  realized  ^ 

An  important  contributor  to  this  dilemma  is  the  still  unexplained  difference 

in  the  crack  propagation  behavior  when  crack  growth  is  initiated  under  impact 

loading  rather  than  under  conventional  quasi-static  conditions.  Specifically, 

(2) 

as  reported  by  Kanninen  et  al,  the  use  of  the  values  obtained  for 
4340  steel  in  quasi-static  initiation  gave  decidedly  poor  agreement  when  used 
to  predict  the  crack  length-time  data  obtained  in  an  impact  test.  In  fact, 
the  Kid  value  inferred  from  the  latter  test  was  roughly  a  factor  of  two 
greater!  Added  to  the  geometry-dependence  that  cast  doubt  on  the  validity 
of  Kid  *  K  D  (V)  as  a  unique  material  property,  there  is  some  concern  about 
the  presently  accepted  elastodynamic  treatments  of  fast  fracture.  This 
paper  describes  a  first  step  towards  a  possible  resolution  of  these  difficulties 
via  the  development  of  an  elastoplastic  dynamic  finite  element  model  for 
the  future  treatment  of  fast  fracture  problems. 

Besides  providing  a  more  realistic  model  of  the  specimen  geometry 
and  the  boundary  conditions,  a  finite  element  method  is  particularly  suitable 
fdr  modeling  nonlinear  material  behavior.  To  avoid  the  use  of  an  extremely 
small  mesh  size  in  the  evaluation  of  tho  dynamic  stress  intensity  factor,  a 
conservation  integral,  J  ^can  be  utilised.  The  '^-integral  is  essentially 
a  consequence  of  several  attempts^  ^ aimed  at  extending  the  regime  of  applicability 
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of  the  path-independent  contour  integral  to  include  dynamic,  elastic- 

plastic,  body  force,  and  thermal  contributions  to  the  energy  release  rate 
under  mixed-mode  conditions.  Then,  crack  propagation  can  be  simulated  via  a 
gradual  crack  tip  nodal-force  release  technique  using  either  crack-length  vs 
time  data  (generation-phase  analysis)  or  a  given  fracture  criterion  (application- 
phase  analysis).  Finally,  to  obtain  more  realistic  material  modelling,  a  strain- 
rate  independent  constitutive  relation  based  on  a  Von  Mises  plasticity  potential 
with  an  isotropic  hardening  rule  has  been  included. 

In  this  paper,  a  general  background  discussion  of  current  problems 
in  dynamic  fracture  is  followed  by  a  description  of  the  salient  features 
of  the  finite  element  based  computational  procedure.  The  validity  of  the 
computational  procedure  is  ascertained  by  solving  problems  for  which  reliable 
analytical  or  experimental  results  are  available.  Results  for  both  stationary 
and  propagating  cracks  are  presented.  Also  presented  is  a  comparison  of  the 
results  of  linear  elastodynamic  vs  elastoplastic  dynamic  analyses  performed 
on  a  three-point-bend  specimen  of  AISI  4340  steel  under  impact  loading 
conditions. 


BASIS  OF  THE  COMPUTATIONAL  PROCEDURE 


Background 


Until  recently  there  has  been  a  controversy  between  the  use  of 
static  or  dynamic  treatments  for  the  arrest  of  a  rapidly  propagating  crack. 

On  the  basis  of  results  obtained  by  Hahn  et  al^  and  Kalthoff^,  it  is  now 
widely  believed  that  a  dynamic  based  analysis  is  the  more  correct.  Nevertheless, 
workers  in  dynamic  fracture  mechanics  are  now  faced  with  other  problems. 
Analytical  studies,  which  have  been  until  recently  based  predominantly  on  linear 
elastodynamic  analyses,  have  brought  increased  understanding  of  the  problem. 

But,  they  have  also  brought  to  light  some  new  problems. 
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(2) 

The  problem  identified  by  Kennlnen  et  *1  concerns  the  marked 
differences  in  the  initiation  and  growth  of  cracks  initiated  under  different 
loading  conditions.  Their  key  result  is  shown  in  Figure  1.  This  experiment 
was  performed  on  an  impact  loaded  three-point-bend  specimen  of  AISI  4340 
steel  (yield  strength  a  200  ksi) .  It  can  be  seen  that,  by  using  values 
obtained  from  quasi-static  initiation  tests  (i.e.,  ■  65  +  .044  V),  a 

very  poor  prediction  is  obtained.  Instead,  the  value  ■  170,  which  has  no 
apparent  connection  with  the  "established"  value,  is  needed  for  good  agreement. 

Because  the  analytical  results  in  this  study  were  obtained  by  a 
relatively  simple  elastodynaalc  finite  difference  scheme,  one  possible  reason 
for  the  discrepancy  would  be  the  analysis  procedure  itself.  Solving  the  same 
problem  with  an  entirely  different  -  and  preferably  improved  -  analytical 
procedure  should  remove  any  such  doubt.  For  this  purpose  a  finite  element 
computer  code  with  isoparametric  element  formulation  was  developed.  To  further 
depart  from  the  previously  used  global  energy  balance  approach  for  the 
calculation  of  stress  intensity  factor,  the  J-integral  approach  was  implemented. 

A  second  possible  reason  for  the  discrepancy  could  ba  the  assumption 
of  linear  elastic  material  behavior.  The  finite  element  code  was,  therefore, 
enriched  to  model  the  material  behavior  in  accordance  with  a  user  supplied 
uniaxial  stress-strain  curve.  This,  and  other  reasons  for  the  discrepancy, 
will  be  dealt  with  in  a  later  section. 

Outline  of  the  Solution  Procedure 

The  approach  used  for  the  solution  of  the  equations  of  motion 
in  the  analytical  procedure  presented  here  employs  a  displacement-based 
finite  element  method.  An  isoparametric  finite  element  formulation  with 
linear  and  quadratic  shape  functions  in  a  two-dimensional  space  is  used. 


Crock  Growth,  mm 
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FIGURE  1.  CRACK  GROWTH  IN  4340  STEEL  UNDER  IMPACT  LOADING  I 

L  -  0.181  in,  a  -.0095  in,  B  -  .0158  m,  W  -  .038  m  j 
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General  quadrilateral  e lament ■  with  e  variable  number  of  node*  in  both 
plane-stress  and  plane-strain  conditions  may  be  used.  If  so  desired  the 
l//r  or  1/r  stress  singularity  at  the  crack  tip  may  be  imposed  by  using 
the  quarter  point  approach^^),  However,  this  feature  was  not  util¬ 
ised  in  any  of  the  analyses  presented  in  this  paper. 

The  modified  Newton-Raphe on  approach^ ^  is  used  for  elastic- 
plastic  analysis.  Von  Mises  yield  condition  with  isotropic  strain  har¬ 
dening  is  assumed.  For  time  integration,  either  an  explicit  (central 
difference)  or  an  implicit  (Newmark-Bata)  scheme  may  be  used.  Due  to  it 
being  Inherently  more  stable,  the  Implicit  scheme  offers  computational 
advantages  in  the  solution  of  dynamic  fracture  problems.  All  the  re¬ 
sults  included  in  the  present  paper  were  obtained  by  using  the 
Newmark-Beta  time  integration  schema^). 

Crack  growth  is  modelled  by  gradually  releasing  the  force  ex¬ 
perienced  by  the  crack  tip  at  a  given  Instant  of  time  in  several  steps. 
Details  of  the  crack  growth  modelling  scheme  are  desw'lbed  in  a  paper  by 
Jung,  Ahmad,  Kanninen,  and  Popaiar^13^.  The  scheme  allows  for  model¬ 
ling  crack  growth  in  both  generation  and  application  phases  of  analysis. 
For  application-phase  analysis  where  the  crack  tip  is  advanced  according 
to  a  selected  fracture  parameter,  a  choice  of  fracture  parameters  is 
necessary.  Currently,  crack-tlp-openmg  displacement  (CTOD),  crack-tip- 
opening  angle  (CTOA),  Mode  I  and  Mods  II  dynamic  stress  intensity  fac¬ 
tors  (Ki  and  Kji) ,  and  the  J  family  of  conservation  integrals^3) 
are  available.  Since  CTOD  and  CTOA  are  obtained  directly  from  the 
finite  element  displacement  solution,  and  Kj  and  Kxi  *r*  obtained  by 

A 

first  calculating  J  for  llnaar  elastic  material,  only  a  description  of 

A 

the  J-integrals  is  included  in  the  following. 

A 

The  J-Integral 

A 

The  mathematical  details  involved  in  the  derivation  of  J  are 
available  in  a  paper  by  Kishlmoto  et  al^3).  Here,  a  general  expres- 

A 

8 ion  for  J  is  taken  in  a  form  which  makes  the  parameter  physically  more 


6 


tractable  and  convenient  to  implement  in  a  computational  scheme > 
This  is  done  by  defining 


(1) 


Where  Che  lower  case  letter  subscripts  stand  for  the  elastic 
(a),  dynamic  (d),  thermal  (t),  plastic  (p),  and  body  force  (b),  contributions 
to  the  integrals  and  K  (-1,2)  indicates  that  each  term  in  Equation  (1) 
is  a  vector. 

Kishimoto  et  al  ^  define '5s' as  follows: 

'Js-  £  Cos  0+^2  Sin  0  ,  (2) 


where  9  is  the  angle  of  crack  extension  measured  anticlockwise  from  the  crac* 
line  (Figure  2).  The  integrals  thru  £  in  Equation  (2)  may  be  expressed 
as  follows:  * 


Jk.-/rtr.(«.V,i5)ir 


kd 


V 


3u. 


ff  c  4^ 

A  ij  a  xk 


dA 


iL 

-ii  3^ 


\  »  Sf*  izrd A  -  S  1/2  a  T  e..  n.dr 

r  +  r  11  K 


dA 


—  Plane  Strain  - 


E 

1-v 


— -Plane  Stress 


Q  -  Direction  of  Crock  Extension 


FIGURE  2.  CRACK  TIP  COORDINATES  USED  IN  THE  DEFINITION  OF  ^ 


where 


■  clastic  strain  energy  density 

■  traction  vector 

■  stress  tensor 

■  plastic  strain  tensor 
“  displacement  vector 

*»  acceleration  vector 

■  temperature  increment 
•  body  force 


In  Figure  2,  A  is  the  area  enclosed  between  contours  r  and  r<nd 
and  'J' is  defined  as  the  area  Aend  approaches  taro.  For  running  cracks,  it 
is  assumed  that  the  so-called  "process-region"^  3hown  in  Aend  in  Figure  2 
remains  constant  in  dimensions  and  moves  with  the  crack  tip. 

While  there  may  be  some  uncertainty  regarding  the  use  of  'j' as  a 
fracture  criterion,  it  is  highly  appealing  from  a  computational  viewpoint. 
The  strength  of  the  idea  is  in  the  fact  that  other  proposed  forma  of  energy 
release  rate  expressions  such  as  J  of  Rice  J*  of  Blackburn  of 

Eftis  and  Hebowiti^^,  and  expressions  proposed  by  Wilson  and  Yu^\ 


Freund  and  Bui  ^,all  can  be  shown  to  be  specialized  versions  of 

the  ^-integral ^ .  Therefore,  3s  is  at  least  equally  valid  as  a  fracture 

criterion  as  any  of  the  above-mentioned  parameters. 

For  a  linear  elastic  material  (J^  »  0)  and  in  the  absence  of 

body  forces  (J.  ■  0),it  can  be  shown^^  tiUt: 

*b 


(t)  +  Kjj(t) 


1  +  -i_  k2 

J  2y  K1 


%  •  fr 


where  u  is  the  shear  modulus,  and  k  ■  3-4v  for  plane  strain  and 

(3-v)/(l+v)  for  plane  stress. 


(4) 
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In  a  single  mode  situation  the  appropriate  stress  intensity  factor  can 
be  readily  obtained  from  Equation  (3) . 

NUMERICAL  RESULTS 

Solutions  to  some  representative  problems  solved  by  using  the 
analytical  procedure  described  above  are  now  presented.  The  first 
four  problems  were  chosen  primarily  for  ensuring  the  validity  of  the 
solutions  procedure  by  comparing  the  results  with  available  analytical 
solutions  and  with  experimental  results.  In  the  second  and  the  fourth 

(2) 

problems, comparisons  with  the  previously  used  finite  difference  scheme 
are  also  made. 

The  fifth  and  sixth  problems  were  selected  to  demonstrate  the 
differences  in  obtained  by  using  quasi- static  analysis  vs  the  dynamic 
analysis,  and  to  illustrate  the  effect  of  plasticity  even  in  a  high 
strength  material,  AISI  4340  steel.  Note  that,  in  all  cases  involving  linear- 
elastic-material  assumption,  the  stress  intensity  factors  were  calculated 
by  the  'j'  approach. 

Problem  1:  Stationary  Crack  in  an  Impulsively-Loaded  Center-Cracked  Panel 

The.  problem  considered  here  is  of  a  center-cracked  plate  (Figure  3) 

loaded  dynamically  in  a  suddenly-applied  uniform  tension  a.  The  material 

was  taken  to  be  linear  elastic  (E  ■  200  GPa,  v  ■  0.3)  having  a  density  of 
3 

5000  Kg/m  .  This  problem  has  been  solved  by  a  number  of  other  authors.  Some 
of  these  results  are  shown  in  Figure  3  along  with  the  results  of  the  present 
analysis.  The  good  correlation  that  is  evident  Indicates  that  the  present 
finite  element  procedure  with  the  ^-integral  provides  sufficiently  accurate 
dynamic  stress  intensity  factors  for  stationary  cracks. 

Problem  2:  Unrestrained  Impact-Loaded  Bend  Specimen 

A  bend  specimen  totally  unrestrained  (no  supports)  is  considered 
(Figure  4) .  In  an  actual  experiment  the  specimen  was  struck  by  a  hammer 
at  an  average  velocity  of  6.88  m/sec  and  was  allowed  to  fly  freely.  In 
Figure  4, the  variation  of  the  dynamic  stress  intensity  factor  with  time. 


Normalized  Dynamic  Siress  Intensity  Factor, 
Kj/la./ira  ) 


Time,  fi  sec 


FIGURE  4.  RESULTS  FOR  UNRESTRAINED  1  SPECIMEN 


L  -  0.181  m,  W  -  0.038  a, 
a  *  0.0095  m,  B  *  0.0158  tn 
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as  obtained  by  the  present  analysis,  Is  compared  with  obtained  by  both 

the  finite  difference  calculation  (2)  and  with  experimental  strain  measurements. 

The  material  properties  used  In  this  analysis  were  those  of  AISI  4340  steel. 

Good  agreement  Is  obtained  which  demonstrates  that  the  'j'-integral  approach 
provides  essentially  the  same  values  for  dynamic  stress  intensity  factors 
as  the  global  energy  balance  approach  used  in  the  finite  difference  scheme 
and.  In  view  of  the  agreement  with  the  strain  gage  data,  that  both  are  correct. 

Problem  3:  Propagating  Crack  in  an  Infinite  Medium 

The  problem  of  a  crack  propagating  at  a  constant  velocity  of  0.2 
times  the  shear  wave  speed  (C  )  in  a  center -cracked  panel  was  analyzed  for 

3 

an  initial  crack  length  to  specimen  width  ratio  (a/w)  of  0.2.  The  results 

/  lg\ 

are  compared  with  analytical  solution  by  Broberg  '  '  in  Figure  5. 

The  results  indicate  that  the  ^-integral  provides  an  effective  means 
of  calculating  dynamic  stress  intensity  factors  for  propagating  cracks. 

However,  it  should  be  noted  that  the  crack  velocity  in  this  comparison  was 
held  constant;  NB,  analytical  solutions  for  a  crack  propagating  at  changing 
velocity  do  not  exist.  In  the  next  problem,  such  a  comparison  is  made  with 
experimental,  as  well  as  numerical  results,  obtained  by  finite  difference 
method. 

Problem  4:  Application -Phase  Analysis  of  a  Propagating  Crack  in  a  Bend  Specimen 

An  elastodynamic  crack  growth  analysis  in  a  three-point-bend 

specimen  of  AISI  4340  steel  was  performed.  Experimental  as  well  as  finite 

difference  analysis  results  for  this  problem  were  first  obtained  by  Kannlnen 
(2) 

et  al  using  K_  ■  65  +  0.044  V  as  a  fracture  criterion.  In  the  present 
tD 

finite  element  computations  the  same  criterion  was  used.  Figure  6  shows  the 
specimen  geometry  used  in  the  analysis  and  a  comparison  of  the  new  results 
with  those  of  Reference  (2). 

This  application-phase  analysis  does  indicate  an  equivalence  between 
tne  ^integral  approach  of  calculating  Kj  and  the  approach  used  in  the  finite 
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FIGURE  6.  A  COMPARISON  OF  FINITE  ELEMENT  ANALYSIS,  FINITE  DIFFERENCE 
ANALYSIS  AND  EXPERIMENT  FOR  A  QUASI -STATICALLY  INITIATED 
CRACK  IN  A  4340  DYNAMIC  TEAR  SPECIMEN. 
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difference  scheme  even  for  e  crack  propagating  at  a  nonconstrnt  speed. 

The  results  also  provide  an  Increased  level  of  confidence  In  the  accuracy 
of  numerical  results. 

Problem  5:  Generation -Phase  Analysis  of  a  Propagating  Crack  in  a  Bend  Specimen 

A  generation  phase  analysis, in  which  an  experimentally  obtained 
crack  length  vs  time  record  (Figure  7)  from  a  dynamic  tear  test  experiment 
on  HY130  steel  was  provided  as  an  input  to  the  computer  code. was  performed. 

The  elastodynamic  finite  element  analysis  then  gave  values  of  as  a  function 
of  time.  These  are  plotted  in  Figure  8.  Also  shown  are  Kq  values  inferred 
by  using  a  handbook  formula  for  three-polnt-bend  specimen  in  which  inertia 
forces  are  ignored.  As  might  be  expected,  the  dynamic  values  oscillate 
around  a  mean  value  given  by  the  quasi-static  solution. 


Problem  6:  Elastic-Plastic  Dynamic  Analysis  of  a  Stationary  Crack  in  a  Bend 
Specimen 


f 


The  three-polnt-bend  specimen  of  AISI  4340  steel  was  analyzed  under 
an  impact  load  with  and  without  the  elastic  material  assumption.  For  the 
elastic-plastic  dynamic  analysis, the  material  behavior  was  described  by  an 
experimentally  obtained  static  stress-strain  curve  from  a  uniaxial  tension 
test.  Strain  rate  effects  on  material  property  are  not  included  in  the  analysis. 

Figure  9  shows  the  variation  with  time  of  the  crack-opening  displacement 
(COD),  measured  at  0.68  mm  behind  the  crack-tip  for  both  elastic  and  elastic- 
plastic  analyses.  In  Figure  10,  values  obtained  are  plotted  against  time. 

It  can  be  seen  that,  even  in  a  high-strength  material,  the  effect  of  plasticity 
appears  to  be  significant.  The  effect  of  plasticity  seems  to  damp-out  the 
oscillations  in  COD  (which  is  directly  proportional  to  and  /T  in  the 
linear  elastic  case);  cf,  Figure  8. 


FIGURE  8.  COMPARISON  OF  ELASTIC  FRACTURE  TOUGHNESS  VALUES  DEDUCED  BY  QUASI-STATIC 
AND  DYNAMIC  CALCULATIONS  FOR  A  HY-130  DYNAMIC  TEAR  TEST  EXPERIMENT 
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DISCUSSION  AND  CONCLUSIONS 


Uncertainties  now  exist  concerning  the  adequacy  of  elastodynamic 
solution  procedures  for  unstable  crack  propagation  and  crack  arrest  analyses, 
the  most  obvious  reason  for  the  existence  of  these  uncertainties  is  the 
neglect  of  crack  ':ip  plasticity  in  such  formulations.  To  address  this 
possibility,  an  elastic-plastic  dynamic  finite  element  capability  has  been 
developed.  It  has  been  shown  in  this  paper,  by  comparisons  with  existing 
dynamic  fracture  mechanics  solutions  and  with  experimental  data,  that  this 
capability  is  more  than  adequate  lor  its  Intended  purposes. 

Some  preliminary  calculations  comparing  the  elastic  and  the  elastic- 
plastic  predictions  for  a  stationary  crack  in  an  Impact  loaded  bend  specimen 
were  also  obtained.  These  indicate  that  the  effect  of  crack  tip  plasticity 
is  significant  even  for  a  high-strength  material.  The  present  capability 
will  next  be  extended  to  treat  a  propagating  crack  whereupon  an  even  more 
prominent  effect  is  expected  to  be  revealed.  But,  whether  or  not  taking 
direct  account  of  crack  tip  plasticity  during  rapid  crack  propagation  will 
resolve  the  questions  that  now  exist  on  the  geometry  and  initiation  mode 
dependence  of  remains  to  be  seen. 
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